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ABSTRACT 

We present results from a subset of simulations from the “Evolution and Assembly of GaLax- 
ies and their Environments” (EAGLE) suite in which the formulation of the hydrodynamics 
scheme is varied. We compare simulations that use the same subgrid models without re¬ 
calibration of the parameters but employing the standard GADGET flavour of smoothed parti¬ 
cle hydrodynamics (SPH) instead of the more recent state-of-the-art ANARCHY formulation 
of SPH that was used in the flducial EAGLE runs. We And that the properties of most galaxies, 
including their masses and sizes, are not signiflcantly affected by the details of the hydrody¬ 
namics solver. However, the star formation rates of the most massive objects are affected by 
the lack of phase mixing due to spurious surface tension in the simulation using standard SPH. 
This affects the efficiency with which AGN activity can quench star formation in these galax¬ 
ies and it also leads to differences in the intragroup medium that affect the X-ray emission 
from these objects. The differences that can be attributed to the hydrodynamics solver are, 
however, likely to be less important at lower resolution. We also find that the use of a time 
step limiter is important for achieving the feedback efficiency required to match observations 
of the low-mass end of the galaxy stellar mass function. 

Key words: cosmology: theory, methods: numerical, galaxies: formation, galaxies: clusters: 
intracluster medium 


1 INTRODUCTION 

Cosmological hydrodynamical simulations have started to play a 
major role in the study of galaxy formation. Recent simulations are 
able to cover the large dynamical range required to study the large- 
scale structure dominated by dark matter as well as the centres of 
halos where baryon physics dominates the evolution. Comparisons 
of such simulations with observations show broad agreement and 
help confirm the predictions of the ACDM paradigm (e.g. Vogels- 
berger et al. 2014; Schaye et al. 2015). 

Galaxy formation involves a mixture of complex processes 
and the numerical requirements to simulate all of the relevant scales 
are enormous. A direct consequence of this is the need to model 
some of the unresolved processes with subgrid prescriptions. Other 
processes, taking place on larger scales, can in principle be fol¬ 
lowed accurately by numerical hydrodynamics solvers. The shock¬ 
ing of cold gas penetrating halos and the turbulence generated by 
supernova activity within galaxies are examples of the processes 
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that can, in principle, be treated by the hydrodynamics solver. Con¬ 
versely, the accretion of gas onto black holes and the formation 
and evolution of stars are examples of processes that occur on 
scales that are too small to be simulated jointly with the large- 
scale environment. In practice, however, these two categories of 
processes are interleaved and it is hence difficult to demonstrate 
convergence even of the purely hydrodynamical processes. Practi¬ 
tioners are therefore forced to chose a numerical hydrodynamics 
solver that gives accurate results at the resolution of interest. 

Many numerical techniques (e.g. Adaptive Mesh Refinement, 
particle techniques, moving-mesh techniques, mesh-free tech¬ 
niques) have been developed over the years to solve the equations 
of hydrodynamics, each of them coming in different “fiavours”, 
i.e. coming with slightly different equations, assumptions and lim¬ 
itations. For the processes that can be simulated using standard 
numerical solvers, the main question is how the various parame¬ 
ters that enter these hydrodynamics solvers affect the formation of 
galaxies in the simulations. For example, it has been reported that 
different numerical techniques and choices of parameters affect the 
disruption of a cold gas blob in a low-density hot medium, a case di- 
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rectly relevant the accretion of gas and satellite of galaxies (Frenk 
et al. 1999; Marri & White 2003; Okamoto et al. 2003; O’Shea 
et al. 2005; Agertz et al. 2007; Wadsley et al. 2008; Mitchell et al. 
2009; Keres et al. 2012; Sijacki et al. 2012). In principle, the values 
of these numerical parameters can be set by performing controlled 
numerical experiments for which the solution is known. 

In the case of simulations using Smoothed Particle Hydrody¬ 
namics (SPH) solvers (Lucy 1977; Gingold & Monaghan 1977) 
(see Springel 2010a; Price 2012, for reviews), the free parameters 
relate to the treatment of shocks, artificial viscosity and conduction 
and are related to the way the SPH equations are derived, leading to 
different flavours of the technique. Performing controlled tests such 
as Sedov explosions or Kelvin-Helmoltz instabilities (e.g. Price 
2008; Read et al. 2010; Springel 2010b; Hu et al. 2014; Hopkins 
2015; Beck et al. 2015) enables the simulator to identify well- 
motivated values for the parameters and understand the limitations 
of the formulation. Early flavours of SPH had issues dealing with 
discontinuities in the fiuid. One of many examples of this prob¬ 
lem is the “blob test” of Agertz et al. (2007) which was widely 
used in the literature to demonstrate the failure of SPH. A lot of 
effort has then been spent by the community to improve the situ¬ 
ation and many alternative solutions have been proposed to over¬ 
come the appearance of spurious surface tension that prevents the 
correct mixing of phases. Solutions using either an alternative for¬ 
mulation of the equations in which the discontinuities are smoothed 
were proposed (e.g. Ritchie & Thomas 2001 ; Read et al. 2010; Abel 
2011; Saitoh & Makino 2013) as well as solution involving addi¬ 
tional terms diffusing material across the discontinuities (e.g. Price 
2008). Both types of solutions to the discontinuity problem present 
shortcomings (see discussion at the end of Section 2.2) and this mo¬ 
tivated the implementation of SPH used in this paper, which uses a 
combination of both solutions. Although one can in principle cali¬ 
brate the free parameters using tests, it is unclear whether there is 
a single set of values that is suitable for all problems and whether 
these parameter values are also the best choice when performing 
simulations of very hot and diffuse conditions, such as those present 
in the hot haloes of galaxies (e.g. Sembolini et al. 2015). More¬ 
over, the large gap in resolution between these controlled experi¬ 
ments and cosmological simulations makes the extrapolation of the 
solver’s behaviour a difficult and uncertain task. The correct treat¬ 
ment of entropy jumps across shocks or of the spurious viscosity 
that can appear in differentially rotating disks can have direct con¬ 
sequences for the population of simulated galaxies. 

In their comprehensive study of galaxy formation models, 
Scannapieco et al. (2012) used multiple hydrodynamics solvers 
coupled to multiple sets of subgrid models to simulate the forma¬ 
tion of galaxies in a single halo and to study the relative impact 
of the choice of solver and subgrid model. One of their main find¬ 
ings was that the variations in the hydrodynamics solvers led to 
much smaller changes in the final results than did the changes to 
the subgrid model parameters. This was especially the case for the 
prescription of feedback, which can change the final galaxy tremen¬ 
dously (e.g. Schaye et al. 2010; Haas et al. 2013; Vogelsberger et al. 
2013; Crain et al. 2015). A more controlled experiment was per¬ 
formed by Keres et al. (2012), who compared two hydrodynamics 
solvers but used only a simplified model of galaxy formation and, 
apart for the most massive galaxies, found very little difference in 
the galaxy population despite the large gap in accuracy between 
the hydro solvers tested. Their two simulations, however, displayed 
significant differences in the gas properties, especially in the cold 
gas fractions. More realistic subgrid models, especially of feed¬ 
back, are likely to suppress some of these difference. 


Building on those studies, we attempt to quantify the impact of 
the uncertainties in two different implementations of SPH solvers 
on a simulated galaxy population. The EAGLE simulation project 
(Schaye et al. 2015; Crain et al. 2015) uses a state-of-the-art im¬ 
plementation of SPH, called ANARCHY (Dalla Vecchia in prep., 
see also appendix A of Schaye et al. 2015) and the time-step lim¬ 
iter of Durier & Dalla Vecchia (2012). eagle’s subgrid model pa¬ 
rameters were calibrated to reproduce the observed local Universe 
population of galaxies. In this study, we vary the hydrodynamics 
solver. We compare ANARCHY to the older Springel & Hernquist 
(2002) flavour of SPH implemented in the GADGET code (Springel 
2005) and compare the resulting galaxy population to the one in 
the reference EAGLE simulation and to those in simulations with 
weaker/stronger stellar feedback and to runs without AGN feed¬ 
back. Since EAGLE broadly reproduces the observed galaxy popu¬ 
lation, our test is especially relevant and enables us to disentangle 
the effects of the hydro solver from the effects of the subgrid model. 

This paper is structured as follows. In section 2 the EAGLE 
model and the two fiavours of SPH that we consider are described. 
Section 3 discusses the impact of the hydrodynamics solver on the 
simulated galaxies whilst Section 4 presents differences in the gas 
properties of the haloes. A summary of our findings can be found 
in Section 5. 

Throughout this paper, we assume a Planck2013 fiat ACDM 
cosmology (Planck Collaboration et al. 2014) (/i = 0.6777, U6 = 
0.04825, Um = 0.307 and ag = 0.8288) and express all quantities 
without h factors. 


2 THE EAGLE SIMULATIONS 

The EAGLE set consists of a series of cosmological simulations 
with state-of-the-art subgrid models and smoothed particle hydro¬ 
dynamics. The simulations have been calibrated to reproduce the 
observed galaxy stellar mass function, the relation between galaxy 
stellar mass and supermassive black hole mass and galaxy mass- 
size relation at z = 0.1. The simulations also broadly reproduce 
a large variety of other observables such as the Tully-Fisher rela¬ 
tion and specific star formation rates (Schaye et al. 2015), the H 2 
and Hi properties of galaxies (Lagos et al. 2015, Babe et al. sub¬ 
mitted), the evolution of the galaxy stellar mass function (Furlong 
et al. 2015), the column density distribution of intergalactic metals 
(Schaye et al. 2015) and HI (Rahmati et al. 2015) as well as galaxy 
rotation curves (Schaller et al. 2015) and luminosities (Trayford 
et al. 2015). 

The EAGLE simulations discussed in this paper follow 752^ 
4.3 X 10® dark matter particles and the same number of gas parti¬ 
cles in a 50® Mpc® cubic volume from ACDM initial conditions. 
Note that the simulation volumes considered here are a factor of 
eight smaller than the main 100® Mpc® EAGLE run. The mass of 
a dark matter particle is ttidm = 9.7 x 10® M© and the initial 
mass of a gas particle is rug = 1.8 x 10® M©. The gravitational 
softening length is 700 pc (Plummer equivalent) in physical units 
below z = 2.8 and 2.66 kpc (comoving) at higher redshifts. The 
simulations were run with a heavily modified version of the GAD- 
GET-3 A^-body tree-PM and SPH code, last described in Springel 
(2005). The changes include the introduction of the subgrid mod¬ 
els described in the next subsection as well as the implementation 
of the ANARCHY fiavour of SPH, whose impact on the simulation 
outcome is the topic of this paper. In the next subsections we will 
describe the subgrid model used in the EAGLE simulations with a 
special emphasis on those aspects of the model that are directly 
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impacted upon by the hydrodynamic scheme. For the sake of com¬ 
pleteness, we then briefly describe both the standard GADGET and 
ANARCHY flavours of SPH. 


2.1 Subgrid models and halo identification 

Radiative cooling is implemented using element-by-element rates 
(Wiersma et al. 2009a) for the 11 most important metals in the 
presence of the CMB and UV/X-ray backgrounds given by Haardt 
& Madau (2001). To prevent artiflcial fragmentation, the cold and 
dense gas is not allowed to cool to temperatures below those cor¬ 
responding to an equation of state Peos oc that is designed to 
keep the Jeans mass marginally resolved (Schaye & Dalla Vecchia 
2008). Star formation is implemented using a pressure-dependent 
prescription that reproduces the observed Kennicutt-Schmidt star 
formation law (Schaye & Dalla Vecchia 2008) and uses a threshold 
that captures the metallicity dependence of the transition from the 
warm, atomic to the cold, molecular gas phase (Schaye 2004). Star 
particles are treated as single stellar populations with a Chabrier 
(2003) IMF evolving along the tracks provided by Portinari et al. 
(1998). Metals from supemovae and AGB stars are injected into 
the interstellar medium (ISM) following the model of Wiersma 
et al. (2009b) and stellar feedback is implemented by the stochas¬ 
tic injection of thermal energy into the gas as described in Dalla 
Vecchia & Schaye (2012). The amount of energy injected into the 
ISM per feedback event depends on the local gas metallicity and 
density in an attempt to take into account the unresolved struc¬ 
ture of the ISM (Schaye et al. 2015; Crain et al. 2015). Supermas- 
sive black hole seeds are injected in halos above 10 ^°/i“^Mo and 
grow through mergers and accretion of low angular momentum gas 
(Rosas-Guevara et al. 2013; Schaye et al. 2015). AGN feedback 
is performed by injecting thermal energy into the gas directly sur¬ 
rounding the black hole (Booth & Schaye 2009; Dalla Vecchia & 
Schaye 2012). 

The subgrid model was calibrated (by adjusting the intensity 
of stellar feedback and the accretion rate onto black holes) so as to 
reproduce the present-day galaxy stellar mass function and galaxy 
sizes (Schaye et al. 2015). As discussed by Crain et al. (2015), 
this latter requirement is crucial to obtain a galaxy population that 
evolves with redshift in a similar fashion to the observed popula¬ 
tions (Furlong et al. 2015). 

Haloes were identifled using the Friends-of-Friends (FoF) al¬ 
gorithm (Davis et al. 1985) with linking length 0.2 times the mean 
interparticle distance, and bound structures within them were then 
identifled using the SUBFIND code (Springel et al. 2001; Dolag 
et al. 2009). A sphere centred at the minimum of the gravitational 
potential of each subhalo is grown until the mass contained within 
a given radius, i? 2 oo, reaches M 200 = 200 ( 47 rpcr(^)i^ 2 oo/ 3 ), 
where pcv{z) — /8nG is the critical density at the redshift 

of interest. 


2.2 SPH implementations 

All simulations that are compared in this study use modiflcations 
of the Gadget-3 code. We use both the default flavour of SPH doc¬ 
umented in Springel (2005) and the more recent flavour nicknamed 
ANARCHY (Dalla Vecchia (in prep.), see also appendix A of Schaye 
et al. 2015) implemented as a modiflcation to the default code. For 
completeness, we describe both sets of hydrodynamical equations 
in this section without derivations. For comprehensive descriptions 
and motivations, see the review by Price (2012) and the description 


of the alternative formalism by Hopkins (2013, 2015). A formula¬ 
tion of SPH that is similar to ANARCHY is presented in Hu et al. 
(2014). Note that apart from the differences highlighted in this sec¬ 
tion, the codes (and parameters) used for both types of simulations 
are identical. 

2.2.1 Default GADGET-2 SPH 

In its default version, GADGET-2 uses the fully conservative SPH 
equations introduced by Springel & Hernquist (2002). We will la¬ 
bel this “GADGET SPH” in the remainder of this paper and restrict 
our discussion of the model to the 3D case. As in any flavour of 
SPH, the starting point is a the choice of a smoothing function to 
reconstruct held quantities at any point in space from a weighted 
average over the surrounding particles. In the case of gas density, 
at position x^. the equation reads 

Pi = y^mjVK(|xi, (1) 

3 

where W{\r\,h) is the spherically symmetric kernel function. In 
the case of GADGET, the M 4 cubic B-spline function is used and 
reads 


if 0 ^ r ^ I 
if I < r ^ ft 
if r > h. 

The smoothing length hi of a particle is obtained by requiring 
that the weighted number of neighbours 

iVngb = ^7tft? W (|xi - x4, fti) (2) 

3 

of the particle is close to a pre-deflned constant; A^ngb = 48 in 
our case. Note, however, that contrary to what is often written in 
the literature, GADGET deflnes the smoothing length as the cut-off 
radius of the kernel and not as the more physical FWHM of the 
kernel function (Dehnen & Aly 2012). 

The quantity integrated in time alongside the velocities and 
positions of the particles is the entropic function^ Ai = Pi/ 
deflned in terms of the pressure Pi and polytropic index 7 . The 
equations of motion are then given by 




i-6(7'+6(7; 

2 (l-£) 


dvi 

Hr = -X’ 


ViWijihi) + 


^ip: 


Vtjp- 


'3Hj 


( 3 ) 


where accounts for the gradient of the smoothing length. 




dWijihj 

' dh 


( 4 ) 


and Wij {hi) = IV(xi — Xj, ). In the absence of radiative cooling 
or thermal diffusion terms, the entropic function of each particle is 
a constant in time. Only radiative cooling, feedback events (see the 
previous section) and shocks will change the entropic function. 

In order to capture shocks, artiflcial viscosity is implemented 
by adding a term to the equations of motion (eq. 3) to evolve the 
entropic function accordingly: 


^ This quantity is not the thermodynamic entropy s but a monotonic func¬ 
tion of it. 
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5 + 


dt 

3 

i 

with Wij = (Wzj(/ii) + Wij{hj)) and the viscous tensor (11^^) 
and shear flow switch fi deflned below. Following Monaghan 
(1997), the viscous tensor, which plays the role of an additional 
pressure in the equations of motion, is deflned in terms of the par¬ 
ticle’s sound speed, ci — yJ^Pijpi, as 


Uij — 


with the dimensionless viscosity parameter set to the commonly 
used value of a = 2 in our simulations. Finally, to prevent the 
application of viscosity in the case of pure shear flows, the switch 
proposed by Balsara (1995) is used: 


Alongside the density, which is computed in the usual way (Eq. 1), 
two additional smoothed quantities are introduced in this formula¬ 
tion of SPH: the weighted density 


pi = 


1 


a; 


1/7 


(|xi - yij\,hi) 


( 8 ) 




velocities becomes 

(c-j Cj SWij ) Wij 

a , 

(5) 


Pi + Pj 


dvi 

• /^n (^3 ~ W) • (Xi - X^)\ 

mm 0, ^^^^^- — 

(6) 


V |xi-x,| J 




and its associated weighted pressure. Pi = AipJ. Despite hav¬ 
ing the same units as the regular density, its weighted counterpart 
should only be understood as an intermediate quantity entering 
other equations and should not be used as the gas density. Using 
these two new quantities, the equation of motion for the particle 


p 

-^^QijViWij{hi)p 


A 


1/7 


A: 


1/7 p. 


Ap P] 


j i Vj Wij ( hj ) 


(9) 


with the terms accounting for the gradients in the smoothing length 
reading 


|V-v,| + |Vx v,| + 10-4c,//i,’ 

with the last term in the denominator added to avoid numerical in¬ 
stabilities. The divergence and curl of the velocity held are com¬ 
puted in the standard SPH way (e.g. Price 2012). 

2.2.2 ANARCHY SPH 

The first change in ANARCHY with respect to GADGET is the choice 
of kernel function. More accurate estimators for both the held 
quantities and their derivatives can be obtained by using Wendland 
(1995) kernels (Dehnen & Aly 2012). ANARCHY uses the C 2 ker¬ 
nel. This kernel function is not affected by the pairing instability, 
which occurs when high values of A^ngb are used with spline ker¬ 
nels. It reads 


VF(r, h) 


JLf (1-0^(1+ 4 ^) if o^r^h 

2nh^ \ 0 A r > h. 


To keep the effective resolution of the simulation similar be¬ 
tween the two flavours of SPH, we use A^ngb = 58 with this ker¬ 
nel. This yields the same kernel FWHM as obtained for the cubic 
kernel^ with A^ngb = 48. Note, however, that the C 2 kernel only 
exhibits better behaviour than the cubic spline kernel when large 
numbers of neighbours (Nngh ^ 100) are used (Dehnen & Aly 
2012). We use the C 2 kernel with A^ngb = 58 to be consistent with 
both the EAGLE resolution and the hydrodynamics studies of Dalla 
Vecchia (in prep.) who used the same kernel but more neighbours. 

The equations of motion used in the ANARCHY flavour of SPH 
are based on the Pressure-Entropy formulation of Hopkins (2013), 
a generalisation of the earlier solutions of Ritchie & Thomas 
(2001), Read et al. (2010), Abel (2011) and Saitoh & Makino 
(2013). The two quantities carried by particles that are integrated 
forward in time are again the velocity and the entropic function. 


^ Expressing our resolution in terms of the local inter-particle separation 
(Price 2012; Dehnen & Aly 2012) gives p = FWHM(V(r, h))/Ax = 
1.235 for both kernels. 


^ij — 1 



/ hi dPp\ 
y 3pi dhi J 



-1 


( 10 ) 


The use of the smoothed quantities Pi and pi in the equa¬ 
tions of motion smooths out the spurious pressure jumps appearing 
at contact discontinuities in older formulations of SPH (Saitoh & 
Makino 2013; Hopkins 2013). 

As in all versions of SPH, artiflcial viscosity has to be added to 
capture shocks. In the ANARCHY formulation of SPH, this is done 
following the method of Cullen & Dehnen (2010). Their scheme 
is the latest iteration of a series of improvements to the standard 
(Monaghan 1997) viscosity term that started with the proposal of 
Morris & Monaghan (1997) to assign individual viscosities ai to 
each particle. Improving on the work of Rosswog et al. (2000), 
Price (2004) and Wetzstein et al. (2009), Cullen & Dehnen (2010) 
proposed a differential equation for ai that is solved alongside the 
equations of motion (Eqn. 9): 


— 2/Usig^i (o^loCji Q^-j) jhi, (11) 

with I = 0.01 and the signal velocity Usig,i introduced below. The 
local viscosity estimator a\oc,i is given by 

b2 q. 

O^loc,i — <Tmax“2 . 72 c ’ (12) 

'^sig,i ' ^i 

where amax = 2 and Si = max (O, — ^ (V • v^)) is the shock de¬ 
tector. After passing through a shock. Si = 0 and hence aioc,z = 0, 
leading to a decrease in ai. We impose ai > amin = 0.05 to fa¬ 
cilitate particle re-ordering. The signal velocity is constructed to 
capture the maximum velocity at which information can be trans¬ 
ferred between particles whilst remaining positive: 

Usig,z = max ( ^ (ci-f-Cj) - min(0, Vij • Xij) ) , (13) 

\z j 

with -kij — (xi — Xj)/|xi — Xj| and = Vj — v^. 

The individual viscosity coefficients ai are then combined to 
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enter the equations of motion in a similar way as in the GADGET 
formulation. Equations 5 and 6 are replaced by: 


Uii = - 


On Oij {Ci Cj 3Wij ) Wi 
Pi + Pj 

(Vj - V») ■ (Xj - Xj) 


2 

min ( 0, 


(14) 

(15) 


Note that contrary to Hu et al. (2014), we do not implement ex¬ 
pensive matrix calculations (Cullen & Dehnen 2010) for the calcu¬ 
lation of the velocity divergence time derivative entering the shock 
detector Si as we found that using the standard SPH expressions 
was sufficient for the accuracy we targeted. 

The last improvement included in the ANARCHY flavour of 
SPH is the use of some entropy diffusion between particles. SPH is 
by construction non-diffusive (e.g. Price 2012) and does, hence, not 
incorporate the thermal conduction that may be required to faith¬ 
fully reproduce the micro-scale mixing of gas phases. We imple¬ 
ment a small level of numerical diffusion following the recipe of 
Monaghan (1997) and Price (2008). We compute the internal en¬ 
ergies from the entropies and these are then used in the equations 
for the diffusion. The use of the pressure-entropy formalism (eq. 9) 
prevents the formation of spurious surface tension at contact dis¬ 
continuities (Hopkins 2013). This small amount of numerical dif¬ 
fusion allows the particles to mix their entropies at the disconti¬ 
nuity and hence create one single phase (Dalla Vecchia in prep.). 
The diffusion is hence used to solve a numerical problem and not 
to introduce a macroscopic conduction. This results in cluster en¬ 
tropy profiles in agreement with the results from grid and moving- 
mesh codes (see the comparison of Sembolini et al. 2015, which 
includes ANARCHY). We compute the rate of change of the con¬ 
duction using the second derivative of the energy. This means that 
large conduction values aaiff are triggered by discontinuities in the 
first derivative of the energy, not by smooth pressure gradients as 
in (self-)gravitating objects. Moreover, aaiff may take some time to 
increase while the smoothing of the discontinuity decreases its rate. 
Finally, our rate is lowered to only a few percent of the computed 
value for the value of the free parameter /3 employed, contrary to al¬ 
most all the implementations in the literature. This largely reduces 
spurious pressure waves. More specifically, the equation describing 
the evolution of the entropy includes a new term. 


dAi d^. 

“dT 


1 


3 


rrij 

pi + pj 


( 


A 

Pi 


— (16) 
pj 


with the diffusion velocity given by Vdmnj — niax(ci + Cj + 
(vi — Vj) • (xi — Xj) / |xi — Xj I , 0) and the diffusion coefficient 
by <Tdiff,G' — + <Tdiff,j). The individual diffusion coef¬ 

ficients are evolved alongside the other thermodynamic variables 
following the differential equation 


ddiff,i — /5 


hiVi [Pi/n - i)Pi) 
VPi/il-'^)Pi 


(17) 


where, as discussed above, we adopt /3 = 0.01. We further impose 
0 < adiff,i < 1, but note that the upper limit is rarely reached, 
even for large discontinuities. 


2.3 Thermal energy injection and time-step limiter 

A crucial aspect of the stellar feedback implementation used in EA¬ 
GLE and described in Dalla Vecchia & Schaye (2012), is the in¬ 


stantaneous injection of large amounts of thermal energy Au in the 
ISM. This injection is performed by raising the temperature of a 
gas particle by AT = 10^'^ K, a value much larger than the aver¬ 
age temperature of the warm ISM. In the GADGET formulation of 
SPH, this is implemented by changing the entropy Ai of a particle. 
In the case of ANARCHY, the situation is more complex since the 
densities themselves are weighted by the entropies, which implies 
that a change in the entropy will affect both quantities entering the 
equations of motion of all the particles in a given neighbourhood. 
Hence, changing the internal entropy of just one single particle will 
not lead to the correct change of energy (across all particles in the 
simulation volume) of the gas. The thermal energy injected in the 
gas will be different (typically lower) from what is expected by a 
simple rise in Ai, leading to a seemingly inefficient feedback event. 

This problem is alleviated in the EAGLE code by the use of 
a series of iterations during which the values of Ai and pi are 
changed until they have converged to values for which the total 
energy injection is close to the imposed value: 


Ai^ri-\-l 


Pi,n-\-l 


(7 - l)(Moid + Am) 

ri^n 

Pi.nA^rl^ - miW{0, hi)All^ + miW{0, 
4 I /7 


This approximation is only valid for reasonable values of Au 
and only leads to the injection of the correct amount of energy if 
the energy is injected into one particle in a given neighbourhood, 
as is the case in most stellar feedback events. This scheme typi¬ 
cally leads to converged values (at better than the 5% level) in one 
or two iterations. When large amounts of energy are injected into 
multiple neighbouring particles, as can happen in some AGN feed¬ 
back events, this approximation is not sufficient to properly con¬ 
serve energy (across all particles in a given kernel neighbourhood). 
To avoid this, we limit the number of particles being heated at the 
same time to 30% of the AGN’s neighbours. If this threshold is 
exceeded, the time step of the BH is decreased and the remaining 
energy is kept for injection at the next time step. Isolated explo¬ 
sion tests have shown that this limit leads to the correct amount of 
energy being distributed. 

As was pointed out by Saitoh & Makino (2009), the conser¬ 
vation of energy in SPH following the injection of large amounts 
of energy requires the reduction of the integration time-step of the 
particles receiving energy as well as those of its direct neighbours. 
This was further refined by Durier & Dalla Vecchia (2012), who 
demonstrated that energy conservation can only be achieved if the 
time-step of the particles is updated according to their new hydro- 
dynamical state. This latter time-step limiter is applied in both the 
GADGET-SPH and anarchy-SPH simulations used in sections 3 
and 4 of this paper. We discuss its influence on galaxy properties in 
subsections 3.1 and 3.2. 


3 GALAXY POPULATION AND EVOLUTION 
THROUGH COSMIC TIME 

As discussed by Schaye et al. (2015) and Crain et al. (2015), the 
subgrid models of stellar and AGN feedback are only an incom¬ 
plete representation of the physical processes taking place in the 
unresolved multiphase ISM. In particular, because radiative losses 
and momentum cancellation associated with feedback from star 
formation and AGN in the multiphase ISM cannot be predicted 
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from first principles, the simulations cannot make ab initio predic¬ 
tions for the stellar and black hole masses. In a fashion similar to 
the semi-analytic models, the subgrid models for feedback in the 
EAGLE simulations have therefore been calibrated to reproduce the 
z = 0.1 galaxy stellar mass function and the relations between 
galaxy size and mass and between the mass of the central super- 
massive black hole and the galaxy. The details of this calibration 
procedure are described in Crain et al. (2015). In this section, we 
will present the basic properties of our simulated galaxy population 
when the hydrodynamic scheme is reverted to the commonly used 
GADGET-SPH formalism. We will specifically focus on the galaxy 
stellar mass function and galaxy sizes before turning towards the 
star formation rates. 

We stress that the model parameters have not been recalibrated 
when switching our hydrodynamics scheme back to GADGET-SPH. 

3.1 The galaxy stellar mass function 

In Fig. 1, we show the galaxy stellar mass function (GSMF) ai z = 
0.1 computed in spherical apertures of 30 kpc around the centre 
of potential of the haloes. As discussed by Schaye et al. (2015), 
this choice of aperture gives a simple way to distinguish the galaxy 
and the ICL. The blue and red lines correspond to our simulations 
with the ANARCHY and GADGET fiavours of SPH, respectively. We 
use dashed lines when fewer than 10 objects populate a (0.2 dex) 
stellar mass bin and dotted lines when the galaxy mass drops below 
our resolution limit (for resolution considerations, see Schaye et al. 
2015). The two hydrodynamic schemes lead to very similar GSMFs 
with significant differences only appearing at M* > 2 x 10^^ M©, 
where the small number of objects in the volume prevents a strong 
interpretation of the deviation, based solely on that diagnostic. The 
white circles and grey squares correspond to the observationally 
inferred GSMFs from the GAMA (Li & White 2009) and SDSS 
(Baldry et al. 2012) surveys, respectively. The two simulated galaxy 
populations undershoot the break of the stellar mass function by a 
similar amount and are in a similarly good agreement (<0.2 dex) 
with the data. The choice of hydrodynamic solver seems to only 
impact the mass and abundance of the most massive galaxies in 
our cosmological simulations. We re-iterate that there has been no 
recalibration of the subgrid parameters between the GADGET and 
ANARCHY simulations. 

In order to compare the contribution of hydrodynamics uncer¬ 
tainties to the uncertainties arising from the subgrid models, we 
show using green and yellow lines two additional models using 
the ANARCHY flavour of SPH but with feedback from star forma¬ 
tion injecting half and twice as much energy, respectively. These 
simulations are the models WeakFB and StrongFB introduced by 
Crain et al. (2015) and reduced or increased the number of feed¬ 
back events taking place, whilst keeping the amount of energy 
injected per event constant. They have been run in smaller vol¬ 
umes (25^ Mpc^), leading to poorer statistics at the high-mass end. 
These changes in the amount of energy injected in the ISM lead to 
much larger differences in the GSMF than changing the flavour of 
SPH used for the simulation. 

The large impact of variations of the subgrid model for stel¬ 
lar feedback on the simulated population and on single galaxies 
can also be appreciated from the large range of outcomes of the 
different models in the OWLS suite (Schaye et al. 2010; Haas 
et al. 2013) and Aquila projects (Scannapieco et al. 2012). Our 
work, however, uses a higher resolution than was accessible in the 
OWLS suite for 2 : = 0 and contrary to Aquila uses a cosmolog¬ 
ical volume and can hence study the effect of the hydrodynamics 



Figure 1. The 2; = 0.1 GSMF of the L050N0752 simulations using an¬ 
archy SPH (blue line, the EAGLE default) and gadget SPH (red line). 
Curves are drawn with dotted lines where galaxies are comprised of fewer 
than 100 star particles, and dashed lines where the GSMF is sampled by 
fewer than 10 galaxies per 0.2 dex mass bin. Data points show measure¬ 
ments with 1(7 error bars from the SDSS (Li & White 2009, filled squares), 
and GAMA (Baldry et al. 2012, open circles) surveys. The yellow and green 
lines show the GSMF of the L025N0376 simulations with twice weaker 
and twice stronger feedback from star formation respectively, in a smaller 
25^ Mpc^ volume. The differences due to the choice of hydrodynamics 
scheme are smaller than the differences due to uncertainties in the sub-grid 
modelling. 

scheme from dwarf galaxies to group-sized haloes. The study of 
Keres et al. (2012), which compared the AREPO (Springel 2010b) 
and GADGET-SPH hydro solvers but using simple subgrid models, 
came to the same conclusion: the choice of hydrodynamics scheme 
has little impact on the stellar mass function of simulated galaxies 
at intermediate mass, only the most massive objects are affected. 
Interestingly, the differences they observed in high mass galaxies 
are exactly opposite to our findings: the more accurate solver (in 
their case AREPO) produces more massive galaxies than the sim¬ 
ulation using GADGET-SPH. This confirms that the source terms 
arising from the physical modelling of the unresolved processes 
in the ISM, especially the modelling of AGN feedback (see the 
discussion below in Section 4.2), clearly dominate the uncertainty 
budget. 

We now turn to the impact of the time-step limiter on the sim¬ 
ulated galaxy population. As was demonstrated by Durier & Dalla 
Vecchia (2012), the absence of a time-step limiter leads to the non¬ 
conservation of energy during feedback events. The energy of the 
system after the injection is larger than expected. This implies that 
a simulation without time-step limiter will have a spuriously high 
feedback efficiency. In order to test this, we ran a simulation in a 
25^ Mpc^ volume using the Ref subgrid model and the ANARCHY- 
SPH scheme but with the Durier & Dalla Vecchia (2012) time-step 
limiter switched off. Since this simulation volume is too small to be 
representative, it is more informative to study the relation between 
halo mass and stellar mass. 
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Figure 2. The median ratio of the stellar and halo mass of central galax¬ 
ies, as a function of halo mass M 200 and normalised by the cosmic baryon 
fraction at 2 : = 0.1 for both the L050N0752 ANARCHY SPH (blue line) 
and GADGET SPH (red line) simulations. Curves are drawn with dashed 
lines where the GSMF is sampled by fewer than 10 galaxies per bin. The 
1(7 scatter about the median of the anarchy run is denoted by the blue 
shaded region. The solid and dashed grey lines show the multi-epoch abun¬ 
dance matching results of Behroozi et al. (2013) and Moster et al. (2013) 
respectively. The yellow and green lines show the GSMF of the L025N0376 
simulations with twice weaker and twice stronger feedback from star for¬ 
mation respectively. The cyan line corresponds to the simulation using the 
ANARCHY formulation of SPH and the reference subgrid model, but with¬ 
out the time-step limiter. The absence of the time-step limiter artificially 
increases the efficiency of the feedback and has a greater impact than the 
choice of hydro solver. 


In Fig. 2, we therefore show the relation between halo mass 
(M 200 ) and galaxy formation efficiency (M*/M 2 oo) for central 
galaxies at z = 0.1. As for all other figures, the blue and red lines 
correspond to the ANARCHY-SPH and GADGET-SPH simulations 
respectively, both using the time-step limiter. We show the simula¬ 
tion using twice stronger and twice weaker feedback with green and 
yellow lines respectively. These are the same simulations that were 
shown in Fig. 1. The stronger feedback from star formation leads to 
a lower stellar mass formed in a given halo than in the Ref model, as 
was discussed by Crain et al. (2015). As expected from the GSMF, 
galaxy formation efficiency is strongly moderated by the feedback 
parameters. Finally, we show in cyan the simulation using the Ref 
subgrid model but without the time-step limiter. This simulation 
displays a lower stellar mass in a given halo than its counterpart 
using the limiter. This indicates that the feedback was indeed more 
efficient at quenching star formation in that simulation, as expected 
from the analysis of Durier & Dalla Vecchia (2012). This is a purely 
numerical effect that has to be corrected by the use of small time- 
steps in regions where feedback takes place. The simulation vol¬ 
ume considered for that test is too small to contain a large sample 
of haloes hosting galaxies with significant AGN activity. Durier & 
Dalla Vecchia (2012) argued that the larger the energy jump, the 
larger the violation of energy conservation will be when the time- 


step limiter is not used. As the energy injection in AGN feedback 
events is two orders of magnitude larger than for stellar feedback, 
we expect the masses of galaxies with M* > 3 X 10^°Mq (the 
mass range where AGN feedback starts to be important) to be re¬ 
duced, compared to the Ref model, even more than the galaxies for 
which AGN feedback plays no role. Note that the impact of the 
time-step limiter is much larger than the differences due to the hy¬ 
drodynamics solver, but smaller than the effect of doubling/halving 
the feedback strength. 

3.2 The sizes of galaxies 

Crain et al. (2015) showed that matching the observed GSMF 
does in general not lead to a realistic population of galaxies in 
terms of their mass-size relation and mass build up. Alongside 
galaxy masses, galaxy sizes were therefore considered in the EA¬ 
GLE project during the calibration of the parameters of the sub¬ 
grid model for stellar feedback. Crain et al. (2015) demonstrated 
that numerical limitations tend to make feedback from star forma¬ 
tion less efficient at quenching the galaxies if the feedback occurs 
in dense regions of the ISM. This would lead to galaxies that are 
too compact and with a specific star formation rate at low redshift 
that is lower than observed. As a consequence, they also showed 
that selecting model parameters that lead to galaxies with sizes in 
agreement with observational data was necessary to obtain a realis¬ 
tic population of galaxies across cosmic time. Assessing the depen¬ 
dence of the galaxy sizes on the hydrodynamics scheme is, hence, 
crucial. 

In Fig. 3, we show the sizes of the galaxies in both the ANAR- 
CHY-SPH and GADGET-SPH simulations. The observational data 
sets from Shen et al. (2003) (SDSS, grey line and shading) and 
Baldry et al. (2012) (GAMA, white circles) are shown for compar¬ 
ison. The sizes of the simulated galaxies are computed following 
McCarthy et al. (2012). We fit a Sersic profile to the projected, 
azimuthally-averaged surface density profiles. We then extract the 
half-mass radius of the galaxy, i? 5 o, from this profile when inte¬ 
grated to infinity. To match the observational selection of Shen et al. 
(2003), we select only galaxies that have a Sersic index ris < 2.5. 
We use dashed lines where the (0.2 dex) mass bins contain fewer 
than 10 objects and dotted lines for galaxies that are represented 
by fewer than 600 star particles. The la scatter around the mean in 
the ANARCHY-SPH simulation is shown as the blue shaded region 
for the mass bins that are both well resolved and well sampled. The 
GADGET-SPH simulation presents a similar scatter. 

Both simulations reproduce the observed galaxy size-mass re¬ 
lation. The simulated galaxies lie within 0.1 — 0.2 dex of either of 
the two data sets. As was the case for the GSMF, the galaxy sizes 
are unaffected by the specific details of the hydrodynamics scheme. 
This implies that the two hydro schemes have similar energy losses 
in dense gas regions where feedback takes place. Differences much 
larger than this can be seen when the subgrid model parameters 
are varied, even if one requires the GSMF to match observations 
(Crain et al. 2015). Galaxies with M* > 10^^M© display small, 
but not statistically significant, differences with the objects in the 
GADGET simulation being slightly smaller. This is in agreement 
with the findings of Naab et al. (2007) who, using GADGET-SPH, 
produced massive galaxies too compact compared to observations. 

When considering the galaxy masses, we found that not us¬ 
ing the Durier & Dalla Vecchia (2012) time-step limiter led to an 
increase of the feedback efficiency, although the magnitude of the 
effect was small compared to that of doubling the feedback energy. 
As galaxy sizes were our second diagnostic, we also consider the 
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Figure 3. The sizes, at 2; = 0.1, of disc galaxies in the L050N0752 AN¬ 
ARCHY SPH (blue line) and GADGET SPH (red line) simulations and in 
the ANARCHY SPH model without time-step limiter (cyan line). Size, i^so, 
is dehned as the half-mass radius of a Sersic profile fit to the projected, 
azimuthally-averaged stellar surface density profile of a galaxy, and those 
with Sersic index ris < 2.5 are considered disc galaxies. Curves show the 
binned median sizes, and are drawn with dotted lines below a mass scale 
of 600 star particles, and using a dashed line style where sampled by fewer 
than 10 galaxies per 0.2 dex mass bin. The Icr scatter about the median of 
the ANARCHY run is denoted by the blue shaded region. The solid grey line 
and the grey shading show the median and Icr scatter of sizes for Ug < 2.5 
galaxies inferred from SDSS data by Shen et al. (2003), whilst white circles 
with error bars show sizes of blue galaxies inferred by Baldry et al. (2012) 
from GAMA data. All simulations reproduce the 2; = 0.1 galaxy sizes. 

effect of switching off this limiter on the sizes of our simulated 
galaxies. This model is shown as the cyan line in Fig. 3. The oscil¬ 
lations seen in the curve are due to the smaller volume used for this 
simulation. The sizes of the galaxies are very close to or slightly 
larger than the ones in the default simulation. 

Crain et al. (2015) also showed that using more efficient stellar 
feedback leads (among other things) to higher SSFRs, lower pas¬ 
sive fractions and lower metallicities. We have verified that turning 
off the time-step limiter has the same qualitative effects, although 
the differences are small. We will not consider the effect of turning 
off the limiter further in the rest of this paper. 

3.3 The star formation rates of galaxies 

We now turn to the star formation rates of galaxies. This quantity 
was not used in the parameter calibration process of the ANAR- 
CHY-SPH run (i.e. the default EAGLE model) and is an important 
independent diagnostic of the success of the simulation. Further¬ 
more, since the ISM dictates the star formation rates of galaxies, 
changes in the way the equations of hydrodynamics are solved may 
lead to changes in the SFRs. 

In Fig. 4, we show the average star formation rate per unit vol¬ 
ume. The blue and red lines again correspond to the ANARCHY and 
GADGET flavours of SPH, respectively. Observational data from 



Figure 4. The evolution of the cosmic star formation rate density in both 
the L050N0752 ANARCHY SPH (blue line) and GADGET -2 SPH (red line) 
simulations. The data points correspond to observations from Karim et al. 
(2011) (radio), Rodighiero et al. (2010) (24 /xm), Cucciati et al. (2012) 
(FUV) and Bouwens et al. (2012) (UV). The decline in the star formation 
rate density from 2; = 2 to 2; = 0 is less pronounced in the GADGET 
run, leading to a 65 % higher star formation rate density at 2; = 0. For 
comparison, a model without AGN feedback (yellow line) is shown. The 
star formation rate density in that model has a low-redshift slope similar to 
that of the gadget simulation. 

Rodighiero et al. (2010), Karim et al. (201 1), Cucciati et al. (2012) 
and Bouwens et al. (2012) is also shown. Where applicable, the 
data has been corrected for our adopted cosmology and IMF as de¬ 
scribed in Furlong et al. (2015). In agreement with the data, both 
simulations display a rise in the star formation rate density at high 
redshifts and a fall at z < 2. As was discussed by Furlong et al. 
(2015), the constant offset in star formation rate of 0.2 dex 
between the simulations and observations leads to 20% less stars 
being formed over the cosmic history, consistent with the z = 0.1 
GSMF (Fig. 1), whose “knee” the simulations slightly undershoots. 

The simulation using the GADGET version of SPH predicts a 
higher cosmic star formation rate density than its ANARCHY coun¬ 
terpart between redshifts 2 and 6 but this does not lead to a large 
difference in stellar mass formed by z = 2. However, the higher 
star formation rate seen at z < 1 is important and the smaller de¬ 
crease between z = 1 and z = 0 implies a star formation rate 
that is 65% higher by z = 0 in the simulation using the GADGET 
formulation of SPH. This higher star formation rate can be tenta¬ 
tively related to the larger number of high-mass galaxies seen in 
the GSMF of this simulation and could, hence, indicate a lower 
quenching efficiency of the AGN activity in the largest haloes. An 
extreme version of a model with a low quenching efficiency in large 
haloes is given by a model without AGN feedback. Such a model, 
using the ANARCHY flavour of SPH, is shown using the yellow line 
in Fig. 4. The excess star formation at z < 2 is much larger than in 
the GADGET-SPH based run with AGN feedback, but the slope is 
similar and not steep enough compared to the data. 

Whether the excess star formation rate at low redshift is due 
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Figure 5. The median specific star formation rate of star-forming 

galaxies (M* /M* > 0.01 Gyr“^) as a function of stellar mass at 2 ; = 0.1 
in the L050N0752 ANARCHY SPH (blue line) and GADGET -2 SPH (red 
line) simulations. Dashed line styles are used where the simulation is sam¬ 
pled by fewer than 10 galaxies per 0.2 dex mass bin. The Icr scatter about 
the median of the anarchy run is denoted by the blue shaded region. 
Observational data points with error bars correspond to the median and 
1(7 scatter of the SSFR from GAMA by (Bauer et al. 2013, grey circles) 
and SDSS+WISE by (Chang et al. 2015, white squares). Galaxies with 
M* > 2 X 1O^°M0 have a significantly higher specific star formation 
rate in the GADGET SPH simulation than in the ANARCHY SPH one, but 
the decrease is smaller than when AGN activity is turned off (yellow line). 

to large haloes can be confirmed by looking at the specific star 
formation rate (SSFR) of the simulated galaxies. This quantity is 
shown in Fig. 5 as a function of stellar mass. We limit our selec¬ 
tion to star-forming galaxies by excluding objects with M*/M* < 
0.01 Gyr“^. As was the case for the stellar mass of the galax¬ 
ies, we measure the SFR within a 30 kpc spherical aperture. The 
red and blue lines show the mean SSFR in the simulations using 
the GADGET and ANARCHY fiavours of SPH respectively. As for 
other figures, the lines are dashed when a given mass bin is sam¬ 
pled by fewer than 10 objects. The blue shaded region indicates 
the 1-cr scatter in the ANARCHY-based simulation. The GADGET- 
based simulation displays a scatter of the same magnitude. For 
comparison, we show the SSFR inferred from observations in the 
GAMA survey by Bauer et al. (2013) (grey circles) and observa¬ 
tions by Chang et al. (2015) using recalibrated star formation rate 
indicators based on SDSS-iWISE photometry (white squares). Sim¬ 
ulated galaxies with masses ~ 10^^M© are in agreement with 
the Bauer et al. (2013) data, whilst lower-mass objects exhibit a 
specific star formation rate lower than observed with the discrep¬ 
ancy reaching ~ 0.3 dex at M* ~ lO^M©. Schaye et al. (2015) 
showed that part of this discrepancy goes away if the resolution of 
the simulation is increased. Interestingly, the recalibrated star for¬ 
mation tracers of Chang et al. (2015) lead to lower specific star for¬ 
mation rates in excellent agreement with the EAGLE results. Both 
the GADGET and ANARCHY simulations show the same behaviour 
at low masses. 


At the upper end of the mass spectrum the two simulations 
do, however, differ. The star formation rate of galaxies with M* > 
2 X lO^^M© is significantly larger for the GADGET formulation of 
SPH. At M* ~ 10^^M©, the discrepancy is 0.3 dex. 

Complementary to the SSFR of the star forming galaxies, the 
passive fraction provides a good diagnostic of the efficiency with 
which SF is quenched in large galaxies This quantity is shown in 
Fig. 6 for both our simulations. Galaxies are considered passive if 
their SSFR is smaller than 0.01 Gyr“^, which is an order of mag¬ 
nitude below the observed median SSFR for star forming galaxies 
at that redshift. For comparison, the data points show the fractions 
inferred from SDSS data by Gilbank et al. (2010) and Moustakas 
et al. (2013). We only show points for the simulated population 
at masses for which there are at least 100 particles at the median 
SSFR (see Schaye et al. 2015). 

The two simulations present a very different behaviour for 
galaxies with M* > 2 x 10^°M©. Whilst the ANARCHY-SPH sim¬ 
ulation follows the trend seen in the observational data, the GAD- 
GET-SPH simulation shows a constant passive fraction of ~ 15% 
at masses up to M* = 2 x 10^^M©. At larger masses, the frac¬ 
tion is 0, implying that all galaxies are star-forming, in disagree¬ 
ment with the data that indicates that almost all galaxies (> 80%) 
of that mass range are passive. Note, however, that there are only 
20 galaxies with M* > 10^^M© in the simulation volume and 
that the fractions displayed in Fig. 6 are, hence, affected by small 
number statistics. Since the ANARCHY and GADGET simulations 
use the same initial conditions, the comparison between the two 
schemes is, however, still meaningful. Switching from ANARCHY 
to standard GADGET has qualitatively a similar effect as switching 
off AGN feedback (yellow line). 

The shortage of passive galaxies in the GADGET simulation at 
the high-mass end of the galaxy population and the higher SSFR for 
high-mass objects both indicate that the star formation quenching 
processes are inefficient in the largest haloes. This higher star for¬ 
mation rate at low redshift in high-mass haloes leads to an increase 
of the stellar mass of massive galaxies as was hinted at by the dif¬ 
ference in the GSMFs between the two simulations at z = 0.1 (Fig. 
1). AGN feedback, which is the main source of quenching in our 
model for galaxies with M* > 2 • 10^°M©, seems to be insuffi¬ 
ciently effective at quenching star formation in large haloes in the 
GADGET simulation. 

It is worth mentioning that we cannot eliminate the possibil¬ 
ity that a re-calibration of the subgrid parameters could bring the 
GADGET simulation into agreement with the data. By changing the 
frequency of the AGN events or the temperature to which the gas 
is heated during such an event, it might be possible to quench star 
formation in large galaxies even when the GADGET formulation of 
SPH is used. It is, however, unclear if this could be achieved and 
whether subgrid parameters should be used to compensate for the 
shortcomings of a particular hydro scheme. Similarly, simulations 
run at different resolutions might lead to different conclusions (if 
the subgrid parameters are kept fixed). Note that simulations run 
at a lower resolution (such as the low-redshift versions of OWLS 
Schaye et al. (2010) and cosmo-OWLS (Le Brun et al. 2014)) have 
fewer resolution elements in the haloes and may hence not suffer as 
much from the lack of phase mixing (see discussion below). A full 
exploration of the subgrid model parameter space or a comprehen¬ 
sive resolution study are, however, beyond the scope of the present 
paper. 

The effectiveness of the AGN feedback can be related to the 
state of the gas surrounding the galaxies and in the whole halo. The 
difference can be understood as follows. The accretion of cold gas 
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Figure 6. The fraction of passive galaxies < 0.01 Gyr“^) at 

2; = 0.1 in both the L050N0752 ANARCHY SPH (blue line) and GADGET 
SPH (red line) simulations. We only show mass bins that correspond to 100 
or more star-forming particles for the median SSFR. The grey circles and 
white squares correspond to the passive fractions inferred from the SDSS 
data by Gilbank et al. (2010) and by Moustakas et al. (2013). The passive 
fraction is far too low for galaxies with M* > 2 x IO^^Mq in the GADGET 
simulation, in a similar fashion to the ANARCHY simulation without AGN 
feedback (yellow line). 

onto the galaxies from filaments is the key source of fresh material 
from which stars can be formed in those haloes. The AGN will sus¬ 
tain a hot halo in which these filaments will dissolve. It is likely that 
the spurious surface tension that plagues the density-entropy for¬ 
mulation of SPH used in GADGET does not leave the gas in the hot 
halo in a state where the AGN activity can be effective at stopping 
star formation. An example of these issues would be the inability 
for dense gas blobs to dissolve in a hot halo medium (see for in¬ 
stance the “blob test” problem by Agertz et al. 2007), which could 
allow cold pristine gas in filaments to survive the hot bubbles cre¬ 
ated by the AGN activity and feed the galaxy with gas ready to form 
stars. The better phase-mixing ability of the ANARCHY formula¬ 
tion of SPH is more effective at disrupting infalling filaments and 
prevent them from reaching the galaxies, making the AGN-driven 
bubbles effective at stopping star formation. In this scenario, the is¬ 
sue is not that outflows generated by an AGN are unable to sustain 
a hot halo (we will show that hot haloes are present in both cases), 
it is rather the pristine gas that forms clumps that are unstable and 
cool rather than being mixed in. The next section further explores 
the differences in gas properties of the two simulations. 


4 LARGE- AND SMALL-SCALE GAS DISTRIBUTION 

In the previous section, we showed that the masses and sizes of 
galaxies are only marginally affected by the improvements to the 
hydrodynamics scheme made in the ANARCHY flavour of SPH. 
We also showed, however, that the star formation rates of mas¬ 
sive galaxies are significantly affected by these same improvements 


and argued that some of the differences might be directly related to 
the way in which the different SPH schemes treat the gas in large 
haloes. In this section, we explore this possibility by studying the 
state of the gas both outside and inside haloes. We will focus on the 
largest systems, where the dynamical time is similar to or shorter 
than the cooling time of the hot gas, and hence the hydrodynamic 
forces become important. 

4.1 Gas in large-scale structures 

A simple diagnostic of the state of the gas in a simulation is 
the distribution of the SPH particles or grid cells in the density- 
temperature plane. The different components (ISM, IGM, etc.) can 
then be identified and their relative abundance in terms of mass or 
volume estimated. Since the ANARCHY and GADGET formulations 
behave differently when different phases are in contact or in the 
presence of a shock, it is worth analysing the differences created 
by those schemes. In order to minimize the impact of the subgrid 
models on the distribution of the gas, we start by looking at the gas 
in the inter-halo medium, i.e. the gas outside of haloes. Most of the 
gas that is located outside of haloes has had little contact with star 
forming regions or with the winds driven by AGN and star forma¬ 
tion but some of the material might have been enriched early on 
in proto-haloes (e.g. Oppenheimer et al. 2012). We are hence fo¬ 
cusing on the low-metallicity, mostly primordial, gas before it falls 
onto haloes. This should allow us to consider differences driven 
mostly by the two flavours of the hydrodynamics scheme. 

The haloes have been identified using the FoF algorithm and 
are hence typically larger than the commonly given virial radii. This 
ensures that we are not considering particles that are part of any 
resolved haloes. In both our simulations, we only identify haloes 
that have more than 32 particles, effectively imposing a minimum 
halo mass of MfoF = 3.1 x lO^M©. This analysis is resolution 
dependent via the definition of the minimum halo mass resolved 
by the simulation. If the resolution were increased, one would find 
smaller haloes, meaning that some of the particles that we identify 
as being outside of any halo will become part of small haloes. How¬ 
ever, small haloes are unlikely to host large amounts of star forma¬ 
tion and drive enrichment and feedback. As both simulations have 
been run at the same resolution with the same initial conditions, the 
same objects will collapse and form haloes, ensuring that our one- 
to-one comparison is not compromised by the potential presence of 
smaller unresolved structures. 

In Fig. 7, we show the distribution of the gas outside of all FoF 
groups in the density-temperature plane at z = 0 for GADGET-SPH 
(left panel) and ANARCHY-SPH (right panel). The low-density ma¬ 
terial (tth < 10“^ cm“^) is in a very similar state in the two sim¬ 
ulations with an extended distribution of diffuse material spanning 
more than 4 orders in magnitude in temperature. The higher tem¬ 
perature material has been heated by feedback activity and blown 
out of the haloes in both simulations. Differences start to appear at 
intermediate densities (10“^ cm“^ < nn < 10“^ cm“^). A lot 
more mass resides in that regime in the simulation using the GAD¬ 
GET formulation of SPH. Because of the artificial surface tension 
appearing in GADGET-SPH between different phases in contact dis¬ 
continuities, this dense gas is unable to properly mix with the lower 
density, higher temperature material surrounding it. In the ANAR¬ 
CHY simulation, the use of both the Pressure-Entropy formulation 
of the SPH equations and of a (small numerical) diffusion term has 
allowed this dense gas to dissolve into its surroundings. The differ¬ 
ence is even more striking at higher densities (tth > 10“^ cm“^), 
where no gas is present in the ANARCHY simulation, whilst a sig- 
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Figure 7. The mass-weighted distribution of gas outside of collapsed structures in the density-temperature plane. The left panel shows the 2 ; = 0 distribution 
for the GADGET SPH simulation whereas the right panel shows the equivalent distribution for the ANARCHY SPH simulation. The GADGET SPH run displays 
high-density gas on the imposed equation of state whilst there is no gas in the ANARCHY SPH run above a density of nn > 10“^cm“^. Dense star forming 
gas is mixing with the lower-density, higher-temperature medium in the ANARCHY SPH run, whilst the artificial surface tension introduced by the GADGET 
SPH formulation prevents this gas from dissolving and leads to star formation outside of haloes. 


nificant amount is present in the GADGET one. This difference is 
especially important since, depending on its metallicity, some of 
this dense gas may be star-forming. Star formation is hence taking 
place outside of collapsed structures in the simulation using GAD¬ 
GET. Interestingly, this high-density gas also has a high metallicity 
{Z > O.IZ 0 ). This gas has thus been ejected from haloes after 
having been enriched by star formation. In ANARCH Y-SPH, simi¬ 
lar material would likely be dissolved into the surrounding lower- 
density medium, either outside haloes or in winds inside haloes. 

4.2 Extragalactic gas in haloes 

We find that within haloes differences in the density-temperature 
diagram are best quantified by looking at the distribution of star 
forming gas. We define the IntraGroup Medium (IGrM) as the gas 
within i? 2 oo but outside of 30 kpc masks placed at the centre of 
each subhalo. This excludes the gas present in the ISM or close 
to galaxies and should leave us with a reasonable definition of the 
IGrM. 

In Fig. 8 , we show the star formation rate of the IGrM as 
a function of the halo mass M 200 at z = 0.1 for objects ex¬ 
tracted from the ANARCHY simulation (blue squares) and the GAD- 
GET-SPH simulation (red circles). Haloes with masses M 200 ^ 
10^^M 0 have a higher star formation rate in the IGrM in the simu¬ 
lation using the GADGET formulation of SPH than in the ANARCHY 
simulation. The higher fraction of dense gas (tth > 10“^ cm“^) 
in the GADGET simulation leads to a higher IGrM star forma¬ 
tion rate. The specific star formation of the IGrM corresponds to 
5 X 10“^ Gyr“^ in the GADGET simulation and is more than 
an order of magnitude lower (?^ 4 x 10“^ Gyr“^) for ANARCHY. 
Although these values are low when compared to the typical values 
for galaxies (see Fig. 5), the presence of significant star formation 
in the IGrM indicates that the AGN activity or gravitational heating 


is not effective enough at quenching star formation in the largest 
haloes. 

As the haloes in the GADGET-based simulation exhibit more 
star formation in their IGrM, it is interesting to investigate how 
the dense gas is distributed spatially. To this end, we selected the 
most massive halo (M 200 ~ 2 x 1 O^^M 0 ) in both simulations and 
constructed column density maps of the gas. As we are mainly in¬ 
terested in the dense gas and to increase the clarity of the maps, 
we only select gas with nn > 10“^ cm“^. As discussed above, 
the behaviour of the warm diffuse medium is similar for both for¬ 
mulations of the SPH equations and can hence be safely discarded 
here. 

These dense gas column density maps are shown in Fig. 9 for 
the GADGET (left panel) and ANARCHY (right panel) simulations. 
The large dashed circles indicate the position of the spherical over¬ 
density radius, i? 2 oo ~ 1-1 Mpc, whilst the small solid circles 
indicate the innermost 100 kpc, where the effects of the central 
galaxies on the gas will be maximized. We will not consider this 
central region in the remainder of this subsection since, as was dis¬ 
cussed in section 3, in this region the differences due to the hydro 
solver are likely to be smaller than the ones induced by small vari¬ 
ations in the subgrid parameters. 

The difference between the two maps is striking. The halo 
from the GADGET simulation contains a large number of dense 
clumps of gas at all radii, as was found in the simulations of Kauf- 
mann et al. (2009). These clumps can be seen even inside the inner 
100 kpc where feedback from both the AGN and star formation 
might be expected to disrupt them. These nuggets of dense gas also 
accompany the in-falling satellites. The map extracted from the AN¬ 
ARCHY simulation is much smoother and dense gas is found mostly 
in the wakes of in-falling satellite galaxies following their stripping. 
anarchy’s ability to mix phases in contact discontinuity allows 
dense clumps to dissolve into the hot halo, whereas the spurious 
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Figure 8. The star formation rate of the inter-group medium (IGrM), i.e. 
inside the halo but at least 30 kpc from any identified galaxy, as a function 
of halo mass at 2; = 0.1 for the L050N0752 anarchy SPH (blue squares) 
and GADGET SPH (red circles) simulations. The IGrM is forming signifi¬ 
cantly more stars in group-and cluster-mass haloes (M 200 > SxIO^^Mq) 

in the run using the GADGET SPH scheme. 


surface tension that appears between phases in GADGET SPH al¬ 
lows them so survive and perhaps even grow. Since some clumps 
reach densities that exceed the threshold for star formation, some 
of them will increase the SF rate of the IGrM. Here, the flavour of 
SPH has a direct consequence on the observables extracted from 
the simulation. 

Another observable that may be affected by the choice of hy¬ 
drodynamics scheme is the gas fraction. In Fig. 10, we show the 
result of mock X-ray observations of our haloes. Following the 
method described in Le Brun et al. (2014), we realise mock X- 
ray observations of our haloes and, assuming hydrostatic equilib¬ 
rium, infer the halo mass and gas fraction following the same anal¬ 
ysis that is applied to observational data. For comparison, we show 
data from Vikhlinin et al. (2006), Maughan et al. (2008), Sun et al. 
(2009), Pratt et al. (2009) and Lin et al. (2012). We only selected 
clusters at 2 : < 0.25. As was discussed by Schaye et al. (2015), 
the simulation using the ANARCHY flavour of SPH (blue squares), 
the Ref model of EAGLE, overshoots the extrapolated trend seen for 
higher-mass haloes. This indicates either that the amount of X-ray 
gas in these haloes is too high or that the gas is in the wrong ther¬ 
modynamic state. The analysis of a larger simulation volume with 
more haloes overlapping with the observations motivated Schaye 
et al. (2015) to introduce an alternative model (labelled AGNdT9) 
for which the mock-observation inferred gas fractions are in better 
agreement with the trend in the data. This model uses more sparse, 
but also more energetic AGN heating events and is shown in Fig. 
10 using yellow triangles^. 

Interestingly, the EAGLE Ref model using the GADGET version 


of SPH (red circles) yields results that are very similar to the im¬ 
proved AGNdT9 model combined with ANARCHY-SPH. The gas 
fractions are in reasonable agreement with the data. However, the 
analysis of the dense gas maps and the following discussion in¬ 
dicates that this better agreement is mostly accidental and not a 
success of the model. The X-ray inferred gas fractions are driven 
down by a change in the gas mass in the haloes but also by the 
presence of cold and dense gas in the IGrM that does not emit X- 
ray and hence artiflcially reduces the inferred gas masses. The cold 
clumps lead to the star formation seen in Fig. 7. We note, however, 
that these spurious undisrupted clumps of dense gas are unlikely 
to affect simulations of the IGrM done at lower resolution such as 
those of McCarthy et al. (2010) or Le Brun et al. (2014). Spurious 
surface tension, preventing the mixing of phase, only appears when 
O{10) particles are part of a cold gas fragment. In lower-resolution 
simulation such gas blobs are sampled by fewer particles and will 
mix with their environment. 

The signiflcant difference in star formation rates in massive 
haloes seen between the two formulations of SPH can have con¬ 
sequences for quantities that are directly observable. An example 
of such an observable is the /-band luminosity of groups and clus¬ 
ters (e.g. Sanderson et al. 2013). For galaxies with similar masses 
and metallicities (as is the case when comparing matched pairs of 
galaxies extracted from both our simulations), a higher /-band lu¬ 
minosity indicates a younger population of stars and a higher star 
formation rate over the last billion years. In Fig. 11, we show the 
/-band luminosity as a function of halo mass M 5 oo,hse- The val¬ 
ues are computed by generating mock-observations of our haloes 
as described by Le Brun et al. (2014). Their procedure allows us 
to compute the halo mass and radius assuming hydro-static equi¬ 
librium as is done in observations of actual clusters. The (Cousin) 
/-band luminosity is computed within /? 500 ,hse, the overdensity ra¬ 
dius inferred by assuming hydrostatic equilibrium in the analysis 
of the mock observations. For comparison, we show observational 
data taken from Sanderson et al. (2013), Gonzalez et al. (2013) and 
Kravtsov et al. (2014) as well as the SDSS image stacking result 
of Budzynski et al. (2014). In all cases we selected only clusters at 
z < 0.25. 

As expected from the previous analysis of the star formation 
rates, we And that the /-band luminosity in the groups and clusters 
extracted from the simulation using the GADGET flavour of SPH is 
higher than when using ANARCHY. It is also higher than the trend 
extrapolated from observational data as expected from our analy¬ 
sis of the speciflc star formation rates and the passive fractions for 
massive (M* > 10^^M©) galaxies. In the same flgure, we also 
show the group and cluster luminosities extracted from the simu¬ 
lation using the AGNdT9 model and the ANARCHY SPH scheme. 
The /-band luminosity as a function of mass for that model is very 
similar to the one obtained using the Ref model. The differences 
between the GADGET and ANARCHY based simulations are much 
larger. However, as discussed earlier, changing the model parame¬ 
ters for feedback from star formation will have an even larger effect. 


4.3 ISM and CGM gas 

We now turn to the gas inside galaxies or in their direct vicinity. 
The state of this gas will retain some of the properties of the IGrM 
but will also be directly affected by the subgrid models. 


^ We note that the map of the column density of dense gas of the largest anarchy code (Fig. 9, right panel). There is no large pool of dense gas 

halo in this model is very similar to the one using the Ref model and the clumps floating in the halo. 
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Figure 9. Maps of the column density of dense gas (hh > 0.01 cm“®) in the largest haloes (M 200 ss 2 x IO'^^Mq) of the L050N0752 gadget SPH (left 
panel) and anarchy SPH (right panel) simulations. The large dashed circle shows the location of the spherical overdensity radius i^ 200 , whilst the small 
solid circle in the centre encloses the inner 100 kpc. The halo in the GADGET SPH run contains a large number of dense clumps of gas, as was found by 
Kaufmann et al. (2009) in their simulations, while its counterpart in the ANARCHY SPH run displays a much smoother gas distribution. The spurious surface 
tension appearing in the GADGET formulation of SPH makes it difficult for the dense gas stripped from the in-falling satellites to be disrupted and mixed into 
the IGrM. 


We first focus on the cold and dense phase of the gas. With 
the help of careful simulations using radiative transfer, Rahmati 
et al. (2013) showed that cold (T < 10^'^ K) and dense (tth > 
0.01 cm“^) gas is a good proxy for HI gas. They provide a fitting 
function to compute H I, but for the purpose of this paper, setting 
the HI fraction to 1 for all this cold and dense gas is a sufficiently 
good approximation. In Fig. 12, we show the mass function of the 
HI gas in the ANARCHY (blue line) and GADGET (red line) sim¬ 
ulations. We use dashed lines when the mass bins contain fewer 
than 10 objects and dotted lines when the Hi mass corresponds 
to fewer than 300 SPH particles. We measured the HI mass using 
fixed spherical apertures placed at the centre of each subhalo in or¬ 
der to only select the gas in the ISM and circum-galactic medium 
(CGM). As a point of reference, we show the best-fitting Schechter 
functions to the data of Haynes et al. (2011) (ALFALFA survey) 
and Zwaan et al. (2003) (HIPASS survey). 

As expected from the non-disruption of cold gas in the hot 
halo, there is an over-abundance of massive H I objects in the sim¬ 
ulation using the GADGET variant of SPH. Whilst the simulation 
using ANARCHY is in reasonable agreement with the observations, 
the same model using GADGET overshoots the break in the mass 
function and vastly over-predicts the abundance of HI clouds of 
mass Mhi > 10^° M©. Both simulations under-predict the abun¬ 
dance of low-mass (Mhi ^ 2 x 10^ M©) Hi clouds. As is shown 
by Crain et al. (in prep) for ANARCHY, this is a resolution effect. 
Simulations run with both flavours of SPH exhibit the same be¬ 
haviour in that regime and can then likely be rescued in a similar 
way by increasing the resolution. 

The discrepancy at the high-mass end is another sign that the 
densest gas clumps found in the group- and cluster-like haloes are 
not disrupted by the hot halo. They also seem to survive AGN activ¬ 
ity and the effect of stellar feedback. These large pools of cold gas 


in massive haloes are not observed and are likely to be responsible 
for the spurious star formation seen in the largest galaxies (Figs. 
5 and 6). We note that it might be possible to modify the AGN 
subgrid model so as to disrupt those clouds without breaking other 
constraints imposed on the model. However, it seems unlikely that 
this purely numerical issue can be completely alleviated. Further¬ 
more, the abundance of spurious cold clumps will increase with the 
resolution (as larger fluctuations in the density distribution can be 
sampled), implying that the AGN activity needed to suppress them 
would also have to be modified. 


5 SUMMARY & CONCLUSION 

The aim of this study was to investigate the effects of the improved 
hydrodynamics solver and time stepping used for the EAGLE suite 
of cosmological simulations (Schaye et al. 2015; Crain et al. 2015). 
By running the same simulation without re-calibrating the subgrid 
model parameters with both eagle’s anarchy and the standard 
GADGET formulations of the SPH equations, we were able to iso¬ 
late the effects of the hydrodynamics solver. Thanks to the use of 
the pressure-entropy formulation of SPH (Hopkins 2013), a more 
stable kernel function (Dehnen & Aly 2012), a small amount of 
numerical diffusion (Price 2008), an improved viscosity switch 
(Cullen & Dehnen 2010) and the Durier & Dalla Vecchia (2012) 
time-step limiter, the ANARCHY flavour of SPH is able to reproduce 
a large set of hydrodynamical tests more accurately than the GAD¬ 
GET flavour (Dalla Vecchia in prep., Sembolini et al. 2015). Here 
we investigated whether the better mixing of gas phases implied 
by these changes, as well as the improved treatment of viscosity 
in shear flow, have consequences for the simulation of haloes and 
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Figure 10. The z = 0 gas fractions within i? 5 oo,hse as a function of 
M 5 oo,hse inferred from virtual X-ray observations of the L050N0752 AN¬ 
ARCHY SPH (blue squares) and GADGET SPH (red circles) simulations. 
Data points correspond to measurements from Vikhlinin et al. (2006) (trian¬ 
gles), Maughan et al. (2008) (stars), Sun et al. (2009) (diamonds), Pratt et al. 
(2009) (crosses) and Lin et al. (2012) (pentagons). The ANARCHY SPH Ref 
model overpredicts the gas fractions for group-sized objects but this can be 
solved by using the AGNdT9 prescription for AGN feedback (yellow trian¬ 
gles). The haloes of the gadget SPH run are in better agreement with the 
data as a result of their higher fraction of cold gas that artificially reduces 
the X-ray inferred gas fractions. 


Figure 11.1-Band luminosity within i?5oo,hse as a function of M5oo,hse 
at 2 ; = 0 in the L050N0752 ANARCHY SPH (blue squares) and GADGET 
SPH (red circles) simulations. The yellow triangles show the haloes ex¬ 
tracted from the ANARCHY SPH run with an improved AGN model (AG- 
NdT9).Data points correspond to the observations of Sanderson et al. (2013) 
(triangles), Gonzalez et al. (2013) (stars), Kravtsov et al. (2014) (diamonds) 
and the dashed line represents the SDSS image stacking results of Budzyn¬ 
ski et al. (2014). Where necessary, observations were converted to the I- 
band following Le Brun et al. (2014). The GADGET SPH run overestimates 
the I-band luminosity in the group- and cluster-size objects as expected from 
the absence of passive galaxies in that simulation (see Fig. 6). 


galaxies. Our analysis of the differences can be summarized as fol¬ 
lows: 

(i) Except for the most massive objects, the masses and sizes 
of the simulated galaxies are largely unaffected by the choice of 
SPH flavour. Uncertainties in the subgrid parameters lead to much 
larger differences (Figs. 1 and 3). 

(ii) The absence of the Durier & Dalla Vecchia (2012) time-step 
limiter leads to somewhat more efficient feedback, as expected 
from the non-conservation of energy occurring in feedback events 
when the limiter is neglected. For low-mass galaxies its effect is 
larger than that of the choice of hydro solver but small compared 
to the changes in the subgrid models for feedback (Figs. 2 and 3). 
For AGN feedback the time-step limiter might have a similar or 
stronger effect since the energy per feedback event is greater than 
for stellar feedback. 

(iii) The star formation rates of galaxies in small haloes, where 
the cooling time is smaller than the dynamical time, are unaffected 
by the change of hydrodynamics scheme. However, in massive 
haloes the star formation rates are much higher in the simulation 
using GADGET SPH (Figs. 5, 6 and 11). These differences in 
behaviour can be related to the lower quenching power of the 
AGN activity in that simulation. The lack of phase mixing, coming 
from the spurious artiflcial surface tension appearing at contact 
discontinuities, prevents cold dense gas from dissolving into the 


hot halo (Figs. 7 and 9). 

(iv) This cold dense gas then reaches the central galaxies and leads 
to increased star formation (Figs. 5 and 6) in both the central galax¬ 
ies and intragroup medium (Fig. 8). This also leads to a lower hot 
gas fraction in the haloes (Fig. 10) and an oversestimate of the Hi 
mass (Fig. 12). 

Our results indicate that the improved hydrodynamics scheme 
plays a signiflcant role in hot hydrostatic gas haloes, but not for 
lower-mass galaxies. Our results are resolution dependent and it is 
possible that simulations performed at much higher resolution will 
be more sensitive to the accuracy of the hydrodynamics solver. Fi¬ 
nally, we also stress that some of the differences between the sim¬ 
ulations could potentially be cancelled by changing the values of 
some of the subgrid parameters. 
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by Haynes et al. (2011) and HIPASS data by Zwaan et al. (2003), respec¬ 
tively. The simulation using the GADGET SPH formulation strongly overes¬ 
timates the abundance of massive H i clouds. 
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